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Abstract. We consider two statistical regularities that were used to explain Omori's law 
of the aftershock rate decay: the Levy and Inverse Gaussian (IGD) distributions. These 
^ . distributions are thought to describe stress behavior influenced by various random factors: 
post-earthquake stress time history is described by a Brownian motion. Both distributions 
decay to zero for time intervals close to zero. But this feature contradicts the high immediate 
aftershock level according to Omori's law. We propose that these statistical distributions are 
influenced by the power-law stress distribution near the earthquake focal zone and we derive 
new distributions as a mixture of power-law stress with the exponent i/j and Levy as well 
as IGD distributions. Such new distributions describe the resulting inter-earthquake time 
intervals and closely resemble Omori's law. The new Levy distribution has a pure power-law 
form with the exponent —(1 + iJ)/2) and the mixed IGD has two exponents: the same as 
Levy for small time intervals and — (1 + ip) for longer times. For even longer time intervals 
this power-law behavior should be replaced by a uniform seismicity rate corresponding to 
the long-term tectonic deformation. We compute these background rates using our former 
analysis of earthquake size distribution and its connection to plate tectonics. We analyze 
several earthquake catalogs to confirm and illustrate our theoretical results. Finally, we 
discuss how the parameters of random stress dynamics can be determined through a more 
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detailed statistical analysis of earthquake occurrence or by new laboratory experiments. 
Short running title: RANDOM STRESS AND OMORl's LAW 

Key words: Omori's law; Random stress; Aftershocks; Levy distribution; Inverse Gaussian 
distribution; Seismicity as result of plate-tectonic deformation; Selection bias. 

1 Introduction 

Our aim is to investigate the connection between random stress and the Omori law of af- 
tershock occurrence. Several attempts to explain Omori's law have been published. Since 
stress in the Earth's interiors cannot be easily measured or calculated, these studies usually 
consider stress as a scalar stochastic variable, ignoring for a while its tensorial nature. In 
this work we consider the stress as a scalar. 

Here we use the scalar seismic moment M directly, but for easy comparison we convert 
it into an approximate moment magnitude using the relationship 

m w = jj(log 10 M-9.0), (1) 

(Hanks, 1992), where M is measured in Newton m (Nm). Since we are using almost exclu- 
sively the moment magnitude, later we omit the subscript in m w . 

Kagan (1982) considered stress time history as a Brownian motion, i.e., randomly fluc- 
tuating stress acting on the stressed environment of an earthquake focal zone. Here the 
probability density for the time intervals between earthquake events has a Levy distribu- 
tion with a power-law tail having the exponent —3/2. Kagan and Knopoff (1987) added 
a tectonic drift component to the Brownian motion. For such a case, the distribution of 
inter-event times is the Inverse Gaussian distribution (IGD) which, depending on the value 
of the initial stress drop and velocity of tectonic motion, can exhibit occurrence patterns 
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ranging from the Levy distribution to a quasi-periodic occurrence. 

Matthews et al. (2002) also proposed the IGD as a law for inter-earthquake intervals; 
however, they considered only a limiting, quasi-periodic long-term distribution. The authors 
suggest that "... [the Inverse Gaussian] distribution has the following noteworthy properties: 
(1) the probability of immediate [large earthquake] rerupture is zero; (2) the hazard rate 
increases steadily from zero at t = to a finite maximum near the mean recurrence time ..." 
and real earthquake occurrence follows the same pattern. Ellsworth (1995), Ellsworth et al. 
(1999), Nadeau et al. (1995), and Nadeau and Johnson (1998) present several sequences of 
recurrent quasi-periodic earthquakes which, in their opinion, confirm a smaller coefficient of 
variation than the Poisson process for earthquake recurrence intervals. 

There is ample evidence that after almost any earthquake subsequent events follow 
Omori's law. This has been observed for smaller earthquakes (aftershocks) as well as large 
events comparable in size to the original shock or even exceeding it (Kagan and Jackson, 
1999). Omori's law pattern of power-law rate decaying seismicity is observed at decadal and 
centuries long scales (Utsu et al, 1995; Stein and Liu, 2009; Ebel, 2009). 

Moreover, the presumed pattern of earthquake quasi-periodicity depends on the argument 
that stress drops to zero or close to zero in the focal zone of an earthquake. Thus, it would be 
necessary to wait for some time before a critical stress level is reached. Kagan and Jackson 
(1999) presented several examples of large earthquakes which followed after a short time 
interval (often just a few days) within a focal zone of similar large events. Such an inter- 
earthquake interval is clearly insufficient for stress to be replenished by tectonic motion. We 
present additional evidence for this pattern later in this paper. 

Thus, for relatively short time intervals, earthquakes are clustered in time with the coeffi- 
cient of variation significantly greater than that for the Poisson process (Kagan and Jackson, 
1991). Why is it widely believed that large earthquakes are quasi-periodic in time? Evidence 
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mostly comes from specifically selected sequences (Bakun and Lindh, 1985; Ellsworth, 1995) 
or from paleo-seismic investigations (Ellsworth et al, 1999). However, paleo-seismic inves- 
tigations have poor time resolution: two or more events occurring closely in time are likely 
to be identified as one event. Hence, their coefficient of variation estimates should be biased 
towards smaller values: a more periodic pattern. Some paleo-seismic investigations (Marco 
et al, 1996. Rockwell et al, 2000; Ken- Tor et al, 2001; Dawson et al, 2003) also suggest 
that earthquakes follow a long-term clustering pattern at least in regions of slow tectonic 
deformation. 

Additionally, paleo-seismicity studies the Earth's displacement only at the surface. But 
temporal and spatial distributions of earthquakes rupturing the surface at one site substan- 
tially differ from general earthquake statistics (Kagan, 2005). The earthquake size distri- 
bution is also significantly different for site- and area-based statistics. Rigorously reliable 
statistical properties that are relevant for theoretical studies as well as for seismic hazard 
estimates can only be obtained by analyzing instrumental earthquake catalogs. 

Measurements from instrumental earthquake catalogs indicate that short time intervals 
between large earthquakes are much more frequent than even the Poisson model would 
suggest (Kagan and Jackson, 1991; 1999). Moreover, the problem of biased selection is 
also hard to avoid in historical and paleo-seismic data: sequences not displaying suggested 
patterns or with a small number of events are less likely to be considered. 

Recently, repeating microquakes have been studied in the Parkfield area (Nadeau et al, 
1995; Nadeau and Johnson, 1998, and references therein). These event sequences exhibit 
certain properties of characteristic quasi-periodic earthquakes: regularity of occurrence and 
nearly identical waveforms. They show many other interesting features of earthquake occur- 
rence. However, these micro-events are not characteristic earthquakes in a strict sense: they 
do not release almost all the tectonic deformation on a fault segment while real characteristic 
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earthquakes are assumed to do so. 

As with large characteristic earthquakes, attempts using micro-earthquake quasi- 
periodicity for their prediction have not yet succeeded. Zechar, J. D. (private communication, 
2010) attempted to forecast repeating small earthquakes on the San Andreas fault near Park- 
field, based on event regularity. While retrospective tests indicated that the forecasts based 
on recurrence times were better than random, the forward tests were unsuccessful. 

Moreover, these micro-events apparently occur on isolated asperities surrounded by creep- 
ing zones. Therefore, such an asperity exists as a secluded entity and may produce charac- 
teristic quasi-periodic rupture events, similar to fracture in many laboratory experiments. 
But tectonic earthquakes occur in an environment that is never isolated. This may explain 
their ubiquitous clustering in time and space. 

The quasi-periodic model for characteristic earthquakes (see for example, McCann et al, 
1979; Bakun and Lindh, 1985; Nishenko, 1991) has been put to a prospective test using 
new data. Testing the validity of this hypothesis ended in failure (Rong et al., 2003 and 
references therein; Bakun et al, 2005; Jackson and Kagan, 2006). 

How can we investigate earthquake occurrence patterns in various tectonic zones? Our 
statistical studies of seismicity (see for example, Kagan, 1991, 2007; Kagan and Jackson, 
1991; Kagan et al, 2010) are biased because earthquake rate is dominated by places with a 
high tectonic rate. Thus, the likelihood function maximum is mainly determined by these 
earthquakes. However, we should study earthquake patterns in continental areas (active and 
non-active) where the earthquake rate is low, but the vulnerable human population is large 
(see Stein and Liu, 2009; Parsons, 2009). A naive extrapolation of aftershock sequence rates 
would exaggerate seismic hazard. Instead, we need a convincing tool to produce a truer 
estimate. 

Our recent results suggest that earthquake occurrence can be modeled by the Poisson 
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cluster process: the combination of the long-term rate caused by tectonic strain with short- 
term clustering. For large and medium strain rates, the model appears to work well (see 
again Kagan et al, 2010). The model is producing reasonable long- and short-term forecasts 
of earthquake occurrence both for local (Helmstetter et al, 2006) and global (Kagan and 
Jackson, 2010) seismicity. But to extrapolate such a model to small tectonic rates we would 
need additional arguments; we must extract the rates by comparing seismicity patterns with 
tectonic strain maps. Bird et al. (2010) discuss the dependence of the long-term inter- 
earthquake rate on strain. We need to show how the short-term properties of earthquakes 
behave in the regions of varying strain. As shown below, the transition time from Omori's 
law decay to a quasi- stationary rate clearly increases with a decrease in the strain rate. As 
the result of these studies, we propose that the same statistical model can be applied for 
regions of both high and low tectonic strain. 

2 Aftershock temporal distribution, theoretical analy- 
sis 

2.1 Levy distribution 

Kagan (1982) proposed a heuristical model for an earthquake fracture as follows. At the 
moment that any given earthquake rupture stops, the stress near the edge of the rupture is 
lower than the critical breaking stress for extension of the fracture. The subsequent stress 
history near the earthquake rupture tip depends on other fractures in the neighborhood and 
additional random factors. In this case, the time history of the stress might resemble a 
Brownian random walk. The stress-time function is thus given as a solution to the diffusion 
equation. When the stress reaches a critical threshold level, a new earthquake begins. The 
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distribution density of the time intervals of the 'first passage' (Feller, 1966) is the Levy 
distribution 



where D is the diffusion coefficient, t is the time, and a is the threshold or barrier stress: the 
difference between initial breaking stress and that when rupture ceases. This distribution 
as a function of stress is the Rayleigh law; as a function of time it is the Levy distribution 
(cf. Zolotarev, 1986). In this model, the stress is taken to be a scalar, which corresponds to 
the addition of perfectly aligned stress tensors; the interaction of misaligned tensors will be 
discussed in a later work. The Levy distribution (1) was used by Kagan (1982, see Eqs. 8 
and 9) to model the time dependent part of the space-time sequence of micro-earthquakes. 
Here we are address only the time sequence of events. 

Fig. [1] displays several Levy curves for various values of stress drop, a. The curve's 
tail is a power-law similar to that exhibited by Omori's law. However, for small values 
of time the curves decay to zero: depending on the stress drop, Brownian motion takes 
time to reach a critical stress level. This latter feature, though observed in aftershock 
sequences, is most likely caused by effects from coda waves of both the mainshock and 
stronger aftershocks that prevent identifying smaller aftershocks (Kagan, 2004). This opinion 
is supported by aftershock registration in the high-frequency domain (Enescu et ai, 2009) 
where their time delay is substantially reduced. Measurements of total seismic moment rate 
decay in aftershock sequences (Kagan and Houston, 2005) suggest that aftershocks start 
right after the end (within one minute for a M6.7 earthquake) of a mainshock rupture. 
Therefore, the left-hand decay of the curves in Fig. [T] is not realistic. 

There is another problem in comparison of theoretical distribution with aftershock ob- 
servations. The Levy curves may explain the inter-event time statistics only for the first 
generation (parent-child) of clustered earthquakes. Most observed aftershocks or clustered 




(2) 
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events are generally separated from their 'parents' by multiple generations. For example, 
the present aftershocks of the 1811-12 New Madrid earthquakes or the 1952 Kern County, 
California, earthquake are mostly distant relatives of their progenitors (mainshocks). In the 
Omori law we count all aftershocks, disregarding their parentage, whereas the Levy law de- 
scribes the time distribution for the next event: when the stress level reaches a critical value 
to trigger a new rupture. When we count all events, we decrease the exponent from 1.5 to 
1.0 (Kagan and Knopoff, 1987). 

We assume that the next earthquake may occur in any place in the focal zone of a previous 
event. Attempts to localize the initiation of strong earthquakes have not been successful; for 
example, the 2004 Parkfield earthquake occurred well outside the previously specified area 
(Jackson and Kagan, 2006). 

We also assume that the stress drop after an earthquake has a power-law distribution: it 
is proportional to cr -1- ^, for < if) < 1.0. Kagan (1994), Marsan (2005), Helmstetter et al. 
(2005), Lavallee (2008), Powers and Jordan (2010) support this statement. 

Then for the Levy distribution we obtain 

<f> (t) oc 1 fit) a~ l ~^da = * —77, T ( ^—^-) , (3) 

YKJ J o JKJ 2t^(2tD)^/ 2 \ 2 J ' v; 

where T is a gamma function. The new modified Levy distribution density function is a 
power-law with the exponent —(1 + if)/ 2). The density should be truncated at the left side, 
because aftershocks cannot be observed at short-time intervals close to a mainshock (Kagan, 
2004). Normalizing the distribution would depend on this truncation. Kagan (1991, Eq. 8) 
introduced a minimum time (£m), where M is a seismic moment of a preceding event (see 
also Eq. [T71 below) . Before that time, the aftershock rate is measured with a substantial 
undercount bias (Kagan, 2004). 
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2.2 Inverse Gaussian distribution 

As Matthews et al. (2002) indicate, the Inverse Gaussian Distribution (IGD) has been known 
since 1915. However, this distribution only acquired its name in 1945. This is the name 
found in recent statistical encyclopedias and books (for example, Kotz et al, 2006; Seshadri, 
1999; Chhikara and Folks, 1989), on the Wolfram website, and in Wikipedia, while the 
"Brownian passage-time" (BPT) is little known. Kagan and Knopoff (1987) proposed to use 
this distribution (without calling it IGD) to describe the inter-earthquake time distribution. 
Their major impetus was to explain aftershock statistics (Omori's law) by stress dynamics. 

Here we discuss the appropriate pairwise interval law for a model in which a steadily in- 
creasing source of stress, which we take to be due to plate tectonics, is added to or subtracted 
from a random or diffusion component, where the distribution (J2J) describes the density of 
earthquake recurrence times in the absence of tectonic loading. If the rate of tectonic loading 
is a constant V, the distribution density f(t) is modified to (Feller, 1966, p. 368) become the 
Inverse Gaussian distribution 

(a-Vt) 



m = exp 



2Dt 



(4) 



For tectonic loading velocity V = 0, this equation transforms to fl2]). Fig. [2] displays several 
IGD curves for various values of stress drop. For small values of a and t, the curves are 
similar to the values of Levy distribution (J2J): if the stress drop is insignificant, or the 
tectonic loading influence is small in initial stages of stress behavior. 

Assuming again that the stress in a fault neighborhood is distributed according to a 
power-law, we obtain a new modified distribution of inter-earthquake time, based on the 
IGD (H 



(f>{t) oc / f^a-^da 



o 



exp \-V 2 t/{2D)} 
2t^{2tD)^/ 2 



2' 2D 



+ K,/irfi-^ 1 F 1 (i-t,l,YH 

V D \ 2 1 \ 2' 2' 2D 



(5) 



where \Fi is a Kummer confluent hypergeo metric function (Wolfram, 1999; Abramowitz 
and Stegun, 1972, p. 504). Another expression for the density is 



(6) 



where U a confluent hypergeometric function (ibid.). As in Eq. El both of the above distri- 
butions should be truncated as t — > and can be normalized after it. 

For certain values of ip the expressions ([5] and [6]) can be simplified. For example for if> = 

1 



(f)(t) oc 



2t 



1 + erf V\ 



IJ_ 
2D 



(7) 



For ip = 0.5 we obtain two equations. For a positive V, using Eq. 13.6.3 by Abramowitz and 
Stegun (1972), we transform (jSJ) into 

/ VH 



<p{t) oc 



exp 



1-1/4 



V z t 
AD 



+ h/4 



AD . 



(8) 



2t V AD \ AD 

where J_i/4 and /1/4 are modified Bessel functions (Wolfram, 1999; Abramowitz and Stegun, 
1972, p. 374). 

For a negative V using Eq. 13.6.21 by Abramowitz and Stegun (1972), we transform (jSJ) 
into 

1 / V~ ( T/ 2 A /T/ 2 f\ 

(9) 



0(t) OC 



exp 



K 



2t V vrL) ~ r \ AD J \ AD J ' 
where if 1/4 is a modified Bessel function (Wolfram, 1999; Abramowitz and Stegun, 1972). 

In Fig. [3] we display the new IGD curves for various values of the if) parameter. Although 
the general behavior of the curves remains power-law, the curves change their slope at the 
time value of about 1.0. The curves with V = ±\^D show the distribution difference for 
tectonic loading sign; in the positive case tectonic movement is opposite to the fault displace- 
ment during an earthquake, whereas the negative sign corresponds to motion consistent with 
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the earthquake mechanism. In the latter case, random fluctuations can bring the fault to 
rupture only during the early period of development. 

To show the curves' behavior more clearly in Fig. HJ the PDF values are multiplied by 
£i+i/72_ p or sma u time intervals curves are horizontal, suggesting that the modified IGD is 
similar to the Levy distribution in Eq. EJ 

For small time values, the power-law exponents in Fig. [3] are essentially the same as for 
compounded Levy law ([3]). This can be seen from modified Bessel function approximations 
for small values of the argument (Abramowitz and Stegun, 1972, Eqs. 9.6.7 and 9.6.9), 

/2 n1/4 /2\ 1/4 

i_i /4 (t) oc f-J and K 1/4 (t) oc f-J . (10) 

In the general case the same result can be obtained with Eq. 13.5.5 by Abramowitz and 
Stegun (1972). Then for t — > Eq. [5] transforms into 

(j)(t) oc t- 1 -^ 2 . (11) 

It is obvious from Figs. |3] and H] that, except for the if) = curve, the slope of curves for 
large values of the time increases when compared with t close to zero. Using Eq. 13.1.27 of 
Abramowitz and Stegun (1972) we add an exponential term in ([5]) within the hypergeometric 
function iF\ . Then for t — > oo we obtain 

<f>{t) oct- 1 ^, (12) 

(Abramowitz and Stegun, 1972, Eq. 13.5.1). Therefore, in Eq. |8]the exponent would be —1.5 
(Abramowitz and Stegun, 1972, Eq. 9.7.1), whereas in Eq.[9]the power-law term (Abramowitz 
and Stegun, 1972, Eq. 9.7.2) is multiplied by an exponential decay term 

0(t) oc f-"exp(-^) . (13) 

See Figs. [3] and SJ 
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3 Temporal distribution of aftershocks: Observations 

Many observations of Omori's law behavior have been published (see Utsu et al, 1995 and 
references therein). There are some problems with these measurements. The standard inter- 
pretation of the Omori law is that all aftershocks are caused by a single mainshock. However, 
the mainshock is often followed by strong aftershocks and those are clearly accompanied by 
their own sequence of events, and so on. The second problem is that some earthquakes 
have very few or no aftershocks; such sequences cannot be included in the naive study of 
earthquake sequences, but can be analyzed as stochastic processes. 

Thus, three techniques can be applied to study the temporal distribution of real earth- 
quakes: 1) traditional, phenomenological techniques based on observing individual aftershock 
sequences; 2) using statistical moments of earthquake occurrence, considered as a point pro- 
cess; 3) applying stochastic process modeling to infer the parameter values of earthquake 
temporal interaction. 

3.1 Aftershock sequences 

Beginning with Omori (1894), the temporal distribution of aftershock numbers has been 
studied for more than one hundred years (Utsu et al, 1995). Aftershock rate decay is 
approximated as t~ p with parameter p value close to 1.0 (Kagan, 2004). 

But the simple, superficial study of aftershock rate decay often encounters serious prob- 
lems. First, only relatively rich aftershock sequences can be investigated by direct measure- 
ments; if there are too few aftershocks, their properties can be studied only by combining 
(stacking) many sequences. Second, to isolate individual sequences one should exclude any 
cases when one sequence is influenced by another, an arbitrary procedure which may intro- 
duce a selection bias. Third, an aftershock sequence often contains one or several large events 
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which are clearly accompanied by a secondary aftershock sequence. Taking the influence of 
secondary earthquakes into account is not simple (see Section [3731 for more detail). Fourth, 
some sequences start with a strong foreshock which is sometimes only slightly weaker than 
a mainshock. Again, handling this occurrence presents serious problem. 

Therefore, strong bias may result from directly measuring Omori's law exponents. Two 
other statistical methods considered below enable analysis of the whole earthquake occur- 
rence as a point process, to minimize the problem of data and interpretation technique 
selection bias. 

3.2 Temporal distribution for earthquake pairs 

Kagan and Jackson (1991) investigated space-time pairwise earthquake correlation patterns 
in several earthquake catalogs. They showed that earthquake pairs follow a power-law dis- 
tribution for small time and distance intervals. Kagan and Jackson (1999) showed that 
contrary to the seismic gap model the clustering pattern continues for strong (m > 7.5) 
earthquakes. Large shallow earthquakes can re-occur after a small time interval. 

Table [1] shows the location and focal mechanism difference for m > 7.5 global shallow 
earthquakes in the GCMT catalog (Ekstrom et al, 2005; Ekstrom, 2007) from 1976-2010. 
The table format is similar to Table 1 in Kagan and Jackson (1999). However, here we kept 
only those pairs in the table for which their focal zone overlap (the ^-parameter) is greater 
than 1.0: 

V = , (14) 

where L\ and L2 are the respective rupture lengths for the first and second earthquakes 
in the pair (see Eq. 3 in Kagan and Jackson, 1999) and R is the distance between the 
centroids. Therefore, if 77 > 1.0 the earthquake focal zones would intersect. For several 
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doublets 7] > 2, implying that the smaller event should be largely within the focal zone of 
the larger earthquake. Inspecting the time difference and the 3-D rotation angle between 
focal mechanisms suggests that these high rj pairs may occur after very short time intervals 
and have very similar double-couple mechanisms. 

All earthquakes in the table occur in subduction zones as defined in Kagan et al. (2010). 
However, even with relatively high deformation velocity at these plate boundaries, the inter- 
earthquake time is in most cases substantially lower than the time necessary for tectonic 
motion to restore the critical stress conditions by the occurrence time of the second earth- 
quake (see the last column in Table 1 by Kagan and Jackson, 1999). 

Fig. [5] shows how the normalized number of m > 6.5 shallow earthquake pairs depends 
on the tectonic deformation rate as defined by Bird et al. (2010). Three curves are shown: 
all earthquakes from the GCMT catalog, earthquakes from subduction zones (Kagan et al., 
2010), and events from active continental zones. 

Fig. [6] shows the integration domain for calculating earthquake pair rates. The diagram 
is a square with a side length equal to a catalog duration, T; since the plot is symmetric, 
only the lower-right portion of the square is shown. The first event shown as a filled circle, 
is supposed to be at the square diagonal, the second one at the end of a hatched area. We 
assume that the time difference between earthquakes cannot be less than t (similar to tu 



For the Poisson process, the interval pair density is uniform (see Kagan and Jackson, 
1991, Eq. 1). Thus, the rate is proportional to the hatched area. For the normalized survival 
rate 



where to is the minimum time interval, T is the catalog duration and t is the inter-earthquake 
time interval. 



in Eq.ED. 




(15) 



14 



For the power-law time distribution with distribution density <p(t) oc t 1 e , we obtain 
the normalized survival rate by integrating over the domain shown in Fig. |6j 



In Figs. [TlfTUl the temporal distribution of inter-earthquake times for m > 6.5 event 
pairs is shown as it depends on earthquake rate determined by Bird et al. (2010). Several 
approximations for pair intervals are also given: the Poisson distribution of earthquakes (fl5l) 
in the time span 1976-2010 and a few power-law interval ( fl6|) dependencies. 

Obviously distribution curves consist of two parts: for small time intervals, they follow 
a power-law and for larger intervals the distribution is parallel to the Poisson rate. The 
transition from one behavior to another occurs sooner for zones with a higher tectonic defor- 
mation rate: in Fig. [7] the transition is observed for the time of about 6,000 days; in Fig. [8] 
it is about 4,000-5,000 days, and in Fig. [9] it is less than 3,000 days. 

Fig. HO] shows the time interval distribution for active continental areas. In this case a 
visual inspection suggests that the best approximation is the power-law; no transition to 
the Poisson rate is observable. This absence can be explained by a low deformation rate in 
these zones (see Fig. E]). As was observed for many aftershock sequences in continental and 
slowly deforming areas, aftershock sequences continue according to Omori's law for decades 
and even centuries (Utsu et al., 1995; Stein and Liu, 2009; Ebel, 2009). Here the span of 34 
years covered by the GCMT catalog is likely insufficient to demonstrate the transition from 
an aftershock sequence to a background, Poisson rate. 

We compute an average recurrence time (i) for earthquakes in Figs. ITHTOl The smallest 
value for i is observed for distributions with a significant component of the power-law: Figs. [7J 
and [TOj We also calculate the coefficient of variation (C v ) of earthquake inter-occurrence 
time as a ratio of the standard deviation (a t ) to the average time i (Kagan and Jackson, 




(16) 
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1991). A completely random Poisson occurrence has the coefficient of variation equal to one, 
whereas quasi-periodicity yields a coefficient of less than one. It yields the coefficient larger 
than one for clustered earthquakes. Although the C v estimates are biased when determined 
at a relatively short catalog time span, their mutual relations are indicative of occurrence 
patterns. For Figs. [71TT01 the CVvalues are 1.10, 0.925, 0.841, and 1.977, respectively. These 
evaluations again suggest that earthquakes in areas with a smaller tectonic rate become more 
clustered, and their Poisson component is diminished. 

3.3 Stochastic branching processes 

As we mentioned above, the mainshock is often followed by strong aftershocks and those are 
clearly accompanied by their own sequence of events, and so on. The patterns of multiple 
clustering have been described by stochastic branching processes (Hawkes and Adamopoulos, 
1973; Kagan, 1991; Ogata, 1998). In this model, seismicity is approximated by a Poisson 
cluster process, in which clusters or sequences of earthquakes are statistically independent, 
although individual events in the cluster are triggered. However, Bremaud and Massoulie 
(2001) proposed a model in which all events belong to one branching cluster. 

Kagan et al. (2010) applied the Critical Branching Models (CBM) to analyze earthquake 
occurrences statistically in several global and regional earthquake catalogs. Time intervals 
between earthquakes within a cluster are assumed to be distributed according to a power-law 

^At(At) = 6t e M (At)- 1 - 9 , At>t M . (17) 

This is similar to Omori's law. The parameter 9 is an 'earthquake memory' factor, tu is 
the coda duration time of an earthquake with seismic moment Mj, and /i is the branching 
(productivity) coefficient (Kagan et al, 2010, Eq. 9). 

During a likelihood search the values have been restricted within the interval 0.1 < 6 < 
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1.0; smaller or larger estimates are inadmissible because of physical considerations (Kagan, 
1991; Kagan et al, 2010). Fig. [TT] shows a correlation between two parameters 9 and fi of 
the CBM. Most 0-values are within the interval 0.1 < 9 < 0.5 and are negatively correlated 
with the fi coefficient. 

Several determinations of the time decay exponent are carried out for the ETAS (Epi- 
demic Type Aftershock Sequence) model. Ogata (1998, Tables 2-3) obtained the p-values 
(equivalent to our 1 + 9 coefficient), of the order 1.03-1.14. Ogata et al. (2003, Tables 1-2) 
estimates the p-values to vary within 1.05-1.18. Ogata and Zhuang's (2006, Tables 2-3) 
values are 1.02-1.05. Zhuang et al. (2005) obtained p-values of 1.14-1.15. Helmstetter et al. 
(2006, Table 1) used a branching model similar to the ETAS and calculated p = 1.18 — 1.20. 

To illustrate our fit of the temporal distribution of dependent earthquakes, average num- 
bers of aftershocks following 15 m > 8.0 GCMT earthquakes are displayed in Fig. IT21 (similar 
to Fig. 13 in Kagan, 2004). We use the time period 1977-2003, so that all large earthquakes 
are approximately the same size (8.45 > m > 8.0, i.e., excluding the 2004 Sumatra and 
its aftershocks). Since the GCMT catalog has relatively few dependent events, we selected 
aftershocks from the the U.S. Geological Survey (2010) PDE (Preliminary Determinations 
of Epicenters) catalog. The aftershock rate in the diagram is approximately constant above 
our estimate of the coda duration tu- For the logarithmic intervals, this corresponds to 
the standard form of the Omori law: the aftershock number n a is proportional to 1/t. For 
the smaller time intervals, the aftershock numbers decline when compared to the Omori law 
prediction (1/t). This decline is caused by several factors, the interference of mainshock 
coda waves being the most influential (Kagan, 2004). The decline is faster for weaker events 
(ibid.). 

For theoretical estimates, we used the values of parameters obtained during the likelihood 
function search (Kagan et al, 2010, Table 4) obtained for the full PDE catalog, m t = 5.0: the 
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branching coefficient p, = 0.141, the parent productivity exponent ao = 0.63 (ao = k X 1.5), 
and the time decay exponent = 0.28. Theoretical estimates in Fig. [12]seem to be reasonably 
good at forecasting time intervals on the order of one day. For larger intervals, the expected 
numbers decrease as n a ~ (At) -1 ' 15 : this is stronger than the regular Omori law would 
predict. As we suggested in Section 12.11 the Omori law assumes that all aftershocks are 
direct consequences of a mainshock, whereas a branching model regards any earthquake as 
a possible progenitor of later events. Thus, later aftershocks are the combined offspring of 
a mainshock and all consequent earthquakes. With increase in time, the difference between 
Omori's law and the CBM predictions would increase as well. Marsan and Lengline (2008) 
show that in California catalogs due to cascading aftershock rates for direct and secondary 
triggering differ by a factor of 10 to about 50%. By numerical simulations of the ETAS 
model, Felzer et ai, (2002) and Helmstetter and Sornette (2003) estimated that a substantial 
majority of aftershocks are indirectly triggered by the mainshock. 

The p-parameter in Omori's law is often assumed to be 1.0. This p- value implies that 
the total aftershock number approaches infinity as the duration of the aftershock sequence 
increases. Figs. [3] and HJ however, suggest that in the presence of tectonic loading, the 
time exponent value should increase at longer time periods. Statistical determination of 
the exponent can usually be made for only short periods; thus, we do not yet have a good 
estimate of the p-value for the tail of an aftershock sequence. 
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4 Earthquake statistics and plate-tectonic deformation 



4.1 Stationary earthquake rate due to tectonic deformation 

The tapered Gutenberg-Richter (TGR) relation (Kagan, 2002b) has an exponential taper 
applied to the number of events with a large seismic moment. Its survivor function (1 — 
cumulative distribution) for the scalar seismic moment M is 

F(M) = {Mt/Mf exp ( Mt ~ M ) for M t < M < oo . (18) 



Here M c is the parameter controlling the distribution in the upper ranges of M ('the corner 
moment'), M t is the moment threshold: the smallest moment above which the catalog can 
be considered to be complete; (3 is the index parameter of the distribution; (3 = |6 , b is a 
familiar b- value of the Gutenberg-Richter distribution (G-R, Gutenberg and Richter, 1944). 

By evaluating the first moment of the distribution (ITS]) , we can obtain a theoretical 
estimate of the seismic moment flux (Kagan, 2002b) 

M s = Ml~P T(2 - (3) exp (M /M c ) , (19) 

where T is a gamma function and cto is the seismic activity level (occurrence rate) for 
earthquakes with moment Mo and greater. 

We compare the seismic moment rate with the tectonic moment rate (Mt): 

M t = fiW J J \e\dA = M s / X , (20) 

A 

where \ is the seismic coupling (or seismic efficiency) coefficient, fi is the elastic shear 
modulus, W is the seismogenic width of the lithosphere, e is the average horizontal strain 
rate, and A is the area under consideration. At present, some variables in the equation cannot 
be evaluated with great accuracy; to overcome this difficulty we calculate a product of these 
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variables: the 'effective width' of seismogenic zone W e or coupled seismogenic thickness (Bird 
and Kagan, 2004): 



Bird and Kagan (2004, Eq. 11) propose another, more exact formula for calculating the 
tectonic moment rate appropriate to a plate boundary zone. 

In regions of high seismicity, instead of Eqs. [T9l - [2~i1 we can use measured long-term seismic 
activity to infer the earthquake rate by extrapolating ([18]) to any moment level (see Bird and 
Liu, 2007, Eqs. 4 and 5). Bird et al. (2010) presented an algorithm and tables for a long-term 
world-wide forecast of shallow seismicity based on the Global Strain Rate Map (GSRM) by 
Kreemer et al. (2003). Because GSRM does not estimate tectonic strain-rates of stable plate 
interiors, a simple empirical-averaging method has been used. Thus, the seismicity in plate 
interiors is represented by a uniform rate. 

Since the seismicity level in plate interiors may vary by orders of magnitude, the uniform 
rate may strongly under- or over-estimate the seismicity rate. Therefore, we apply Eqs. IT5H2T1 
to evaluate first the tectonic moment rate and then a long-term forecast for these regions 



By calculating «o for a particular choice of M , we may re- normalize Eq. [TS] and obtain 
earthquake size distribution for any region with a known strain rate and corner moment. 

4.2 Length of aftershock zone 

Eq. [22] above can be used to evaluate the background seismicity level for an aftershock 
zone. Such a calculation would use the area of an earthquake focal zone. This area can be 
estimated by a dimension of aftershock zone for each event. Kagan (2002a) evaluated how the 
aftershock zone size for mainshocks m > 7.0 depends upon on the earthquake magnitude by 



w e = w x - 



(21) 



«o = 



M T xH-P) exp(-M /M c ) 
Ml~ p T(2 - p) 



(22) 
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approximating aftershock epicenter maps through a two-dimensional Gaussian distribution. 
The major ellipse axis is taken as a quantitative measure of the mainshock focal zone size. 

In Fig. [T3] we display the regression curves for GCMT/PDE earthquakes: all earthquakes 
for three choices of focal mechanisms (updated Fig. 6a by Kagan, 2002a). In regression 
curves we use m = 8.25 as a reference point. For example, for the quadratic regression 

L = log 10 I = a + ai (m - 8.25) + a 2 (m - 8.25) 2 , (23) 

where I is the length of the aftershock zone in km. For the linear regression we set a 2 = 0. 
Fig. HH displays the regression in a similar format for active continental tectonic zones (Kagan 
et ai, 2010). 

Table |5] summarizes the results of regression analysis for all global earthquakes, as well 
as for events in subduction zones (trenches) and on active continental zones. Earthquakes 
are also subdivided by their focal mechanism. Other tectonic zones lack a sufficient number 
of m > 7.0 mainshocks to carry out this statistical analysis. 

The following conclusions can be made from Table |2j (a) aftershock zones exhibit similar 
scaling; (b) zone length (£) on average is proportional to moment M 1//3 ; and (c) the value of 
cto parameter (zone length for the m8.25 earthquake) is close to 10 2 ' 5 (316) km for all cases. 
Normal earthquakes (rows 5-6) exhibit slightly different scaling; zone length (£) for the linear 
regression is proportional to moment M 1//2 ' 8 . Scaling for strike-slip earthquakes (rows 7-8) 
also differs a little from average: zone length (£) is proportional to moment M 1 / 3 5 . However, 
the earthquake numbers in these subsets are small, thus it is possible that these variations 
are due to random fluctuations. 

Only three subsets show a substantial nonlinearity: (a) trenches with strike-slip focal 
mechanisms (rows 15-16); (b) continents with all focal mechanisms (rows 17-18); and (c) 
continents with strike-slip focal mechanisms (rows 21-22). However, the earthquake numbers 
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are small in all these plots, and although (b) and (c) aftershocks display zone lengths which 
increase strongly for the largest earthquakes (the feature often quoted in other studies of 
length scaling - see, for instance, Wells, and Coppersmith, 1994, and subsequent publications 
citing that paper), (a) earthquakes exhibit an opposite behavior. 

In all diagrams the standard errors (a) are almost the same for the linear and quadratic 
regression. The maximum errors (e max ) follow the same pattern. This pattern suggest that 
linear regression is sufficient to approximate the data. Although the quadratic regression 
fit yields no statistically significant improvement in almost any diagram, the sign of the 
quadratic correction term is negative for most cases. The negative value of the 02 regres- 
sion coefficient means that increase in the aftershock zone length is weaker for the largest 
earthquakes. This feature contradicts those often quoted in other studies of length scaling 
(see Kagan, 2002a for details). Thus, the slope of the regression curve is either stable or 
decreases at the high magnitude end. No saturation effect for large earthquakes occurs in 
the data. Results in Table [2] imply that the major ellipse axis a (length) of an earthquake 
focal zone can be approximated by 

a = 316 x io( m ~ 8 - 25 )/ 2 km. (24) 

We conclude that earthquake rupture length is proportional to the cube root of moment, 
which implies that width and slip should scale the same way. Otherwise, one of them 
increases less strongly with moment and the other more strongly. For either that would pose 
the problem of "inverse saturation." 

We assume that the majority of aftershocks are concentrated within an ellipse having 2-a 
major axis. The probability that a point lies inside a 2-a ellipse is shown in Eq. 5 by Kagan 
(2002a). If we know the length of an earthquake focal zone, we can calculate its area. We 
assume, for example, that the ratio of the major ellipse axis to the minor axis is 1/4; then 
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area S of the focal zone is 

S = TTO 2 



(25) 



4.3 Example: New Madrid earthquake sequence of 1811-12 

To illustrate the arguments and results of the previous sections, we calculate seismicity 
parameters of the New Madrid earthquake sequence (1811-12) and its consequences. There 
is substantial literature on this sequence (Stein and Liu, 2009; Calais et ai, 2010, and 
references therein). 

Three or four large earthquakes with magnitudes on the order of 7.5-8.0 occurred over a 
few months of 1811-12 in the New Madrid area; aftershocks of these events are still registered. 
As an illustration, we would assume that only one m8 event occurred at that time. If in 
reality earthquakes were smaller than such an event, their total focal zone and combined 
aftershock sequence at the present time would be equivalent to about one m8 mainshock. 

The size of the focal zone can be evaluated by using regression equations in Figs. [13] and 
[HJ The first plot contains many earthquakes but most of these events are in subduction 
zones. The second diagram uses earthquakes in active continental zones; a focal size of 
these earthquakes is likely to resemble the New Madrid area which can be classified as plate- 
interior (Kagan et ai, 2010). Too few large earthquakes are available within plate-interior 
to obtain their features. The difference between regression parameters in Figs. [13] and [14] is 
small; therefore the size of earthquake focal zones either does not change in various tectonic 
zones or changes slightly. 

For an m8 earthquake, calculations yield 227 km and 259 km as the length of the focal 
zone, defined as the 4cr major axis of an ellipse comprising a majority of aftershocks (Kagan, 
2002a). The linear regression is used in both cases: the former value corresponds to Fig. [13] 
and the latter to Fig. [TH The two estimates are similar and roughly correspond to the size 
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of the present aftershock zone, as shown, for example, in Calais et al. (2010). 

To calculate the surface area of an aftershock zone, we assume that the minor axis of 
the ellipse is 1/4 of the major axis, taken as 240 km. Then we obtain the m8 earthquake 
focal £1X621 cLS 11300 km 2 . Taking the average strain rate as e = 10 9 (Calais et a/., 2006), 
we compute the tectonic moment rate ( |20l) : 3.4 x 10 15 Nm/year. Assuming that 50% of the 
tectonic rate is released seismically (Bird and Kagan, 2004), we obtain the background rate 
a = 1.27 x 10~ 3 m > 5 earthquakes per year [we take in fl22|) M c = 10 21 Nm, W = 20 km, 
j3 = 2/3; then T(4/3) = 0.893]. We use ffTB"]) to calculate the recurrence time for an m > 8 
earthquake in the focal zone of the New Madrid events: more than two million years. In 
this computation any m8 earthquake with an epicenter or centroid in the focal zone counts: 
in Eq. (122]) we do not request that the entire rupture of such an earthquake be contained in 
the zone. The recurrence time is an average value; even for events as large as m > 8 the 
earthquake occurrence is clustered (see Section |3~2|) . Thus, a new large earthquake can follow 
after a relatively short time period, as exemplified by the 1811-12 New Madrid sequence. 

We would like to calculate the duration of an aftershock sequence up until the aftershock 
rate decays to the background level. The results in Fig. [12] can be applied to this purpose. 
However, we need to make a correction for the mainshock and aftershock magnitudes (m8 
instead of m8.15 and m > 5.0 instead of > 4.9 in the plot, respectively). In the diagram 
the aftershock rate per one interval (the intervals increase consequently by a factor of two) 
is 7 events. This translates into 4.55 m > 5 events per interval for our choice of magnitudes 
(Kagan at al, 2010). After comparing the background and aftershock rates (we take 4.55 
m > 5 aftershocks per the first day, decaying according to Omori's law, with 1/t rate with 
time), we discover that the aftershock sequence would reach the background rate in about 
3,600 years. This duration estimate agree roughly with Stein and Liu's (2009) value. In 
these calculations, we presume that no independent large earthquake clusters would occur 
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during the aftershock sequence. The possible occurrence of spontaneous events makes any 
evaluation of aftershock sequence duration largely approximate. 

Stein and Liu (2009) obtained aftershock duration values for several sequences using 
Eq. 14 from Dieterich (1994). This equation employs parameters whose values for actual 
earthquake focal zones are not known. Generally, the parameters have been back adjusted 
based on the statistics of earthquake occurrence. This may explain apparently reasonable 
fit of Dieterich's formula to aftershock sequences. 

In contrast, we obtain the aftershock sequence duration by extending the tapered 
Gutenberg-Richter and Omori's laws and using their well-known properties and measured 
geometrical features of tectonic deformation. Moreover, according to Stein and Liu (2009, 
Fig. lc) the New Madrid aftershock rate for the last 50 years was about 0.5 m > 4 events 
per year. Computations based on Omori's law similar to those shown above, yield the rate 
175 years after the mainshock occurrence of about 0.26 m > 4 events per year. This number 
is close to that shown above. 

5 Discussion 

Two classical statistical earthquake distributions largely governed our approach to analyz- 
ing seismicity: Omori's law and the Gutenberg-Richter relation. As explained above, recent 
developments in earthquake size statistics considerably improved our understanding of earth- 
quake occurrence and could lead to significantly better estimates of seismic hazard. For the 
G-R law the earthquake temporal distribution is mostly irrelevant, since size distribution 
of clustered events is largely independent of their history. Similar progress in understand- 
ing earthquake time statistics is much more difficult to achieve. We cannot ignore spatial 
variables and the available data are not as extensive, so the task is more complex. Only by 
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applying rigorous methods, by analyzing carefully systematic and random effects, and by 
critical testing of models and hypotheses we will be able to advance in solving this problem. 

In previous sections we derived the time distribution for earthquake occurrence; the 
distribution is shown to be controlled by power-laws. How can the parameters of these 
distributions be determined? If one excludes the interiors of plates, the tectonic deformation 
rate V is reasonably well known for plate boundaries and for active continents (Kagan et 
ai, 2010; Bird et ai, 2010). The diffusion rate D is presently unknown. If we could obtain 
the earthquake temporal distribution as shown in Figs. [1] and [2], the D evaluation would 
be easy. These distributions are derived for a particular area within an earthquake fault 
zone. However, if the stress in the focal zone of an earthquake is distributed according to 
the power-law with an exponent ip (see Eqs. [3] and [5]), the problem becomes more acute. 

Figs. [3] and @] suggest that the distribution temporal behavior changes when V = \/D. 
This change relates to the first generation of offspring. Thus, we should not be able to see 
it in regular Omori plots which combine many generations of aftershocks. The inversion 
of earthquake occurrence parameters based on stochastic branching processes yields needed 
first generation effects. However, in present models (both CBM and ETAS) as discussed 
in Section 13. 3[ temporal dependence is parameterized by just one exponent. These models 
should in principle demonstrate changes in the temporal pattern, if more complicated tem- 
poral function is applied. However, results from statistical analysis are very uncertain even 
for one-parameter time decay. Given the contemporary quantity and quality of earthquake 
catalogs, it is unlikely more complicated models would be effective in resolving this issue. 

Perhaps new laboratory experiments (Zaiser, 2006) may help solve the problem of diffu- 
sion rate evaluation, but it is not clear whether such measurements are possible. The acoustic 
emission event rate exhibits fore- and aftershock sequences associated with dynamic failure 
of the test specimen (see, for example, Ojala et ai, 2004). These and similar tests can be 
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used to infer the dependence of the Omori law parameters on spatial scale and stress diffusion 
rate. 

Results from statistical analysis of earthquake occurrence in our previous publications 
(Kagan, 2002b; Bird and Kagan, 2004; Kagan et al, 2010), as well as the results reported 
above, suggest that the earthquake process in all tectonic provinces can be described by the 
same model. We advocate the Poisson cluster process with clusters controlled by a critical 
branching process and a power-law time dependence. Combined with the earthquake size 
distribution approximation by the TGR, such a model allows mathematically forecast of a 
spatially variable, time-independent (long-term) earthquake rate. It will optimally smooth 
the seismicity record (Kagan and Jackson, 2010) or translate the plate-tectonic and geodetic 
rate into a seismic hazard estimate (Bird and Liu, 2007; Bird et al, 2010). 

A short-term forecast can be performed by using the temporal properties of earthquake 
clusters, an extrapolation which uses a variant of Omori's law to estimate future earthquake 
rate (Kagan and Jackson, 2010). In Section 14731 we presented an example of such calculations. 

In conclusion, this paper suggests a method for calculating long- and short-term seismicity 
estimates, based on theoretical inference about classical, statistical earthquake distributions: 
the Omori law and the G-R relation. Statistical analysis of earthquake occurrence carried out 
in our previous papers and in this work make such seismic hazard evaluation more reliable 
and accurate. 
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Table 1: Pairs of shallow earthquakes m > 7.5 
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R - centroid distance, $ - 3-D rotation angle between focal mechanisms, At - time interval 
between events, 77 - degree of zone overlap, the ratio of earthquake focal zone sizes to twice 
their distance, see Equations (2,3) in Kagan and Jackson (1999). The total earthquake 
number with magnitude m > 7.50 for the period 1976/1/1-2010/10/25 is 121. The maximum 
epicentroid distance is 250.00 km. 
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Table 2: Aftershock zone log length vs mainshock moment magnitude m 
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ao, a%, <i2 are regression coefficients in Eq. [231 °" is the standard uncertainty; e max is the 
maximum error; n is the number of aftershock sequences. The catalogs time interval is 
1977/1/1-2010/09/21. 
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Fig. 1 




Figure 1: 

Plot of PDFs for the Levy distribution (j2]), D = 1.0; a = 0.001, 0.01, 0.1, 1.0, 10.0, 100.0. 
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Fig. 2 
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Figure 2: 

Plot of PDFs for the IGD distribution @, V = VD; D = 
0.001, 0.01, 0.1, 1.0, 10.0, 100.0, 1000.0. 
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Figure 3: 

Plot of PDFs for the IGD distribution (Eqs. EJ-EJ). 
Curve with squares ([7]): V — a/D, ip = 0.0; 
curve with diamonds OH]): V = VD, ip = 0.5; 
curve with crosses (Q: V = — V^D, ip = 0.5; 
curve with circles (J5J): V = \/D, ip = 0.9. 
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Figure 4: 

Plot of PDFs for the IGD distribution (Eqs. E]-[9j), multiplied by t l+i '' 2 , D = 1.0. 
Curve with squares (JTj): V = a/D, V = 0.0; 
curve with diamonds (jSJ): V = VD, ip = 0-5; 
curve with crosses (J5]): V = — V^D, ip = 0.5; 
curve with circles (J5J): V = \/D, ip = 0.9. 
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Fig. 5 




' All, red (-), pair N = 833, <A> = 39.2, a x = 21 .9 
Trench, green (—), pair N = 752, <A> = 41 .9, a. = 20.9 
Act. cont., blue (-.), pair N = 37, <A> = 12.2, a. =11.7 
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Figure 5: 



Normalized cumulative distribution of earthquake rate for m > 6.5 shallow earthquake pairs 
in all zones, trench (subduction) zones, and active continental zones. The pair numbers N, 
average rate < A > and its standard deviation <j\ are also shown. Earthquake rate is taken 
from a table by Bird et al. (2010) for a magnitude threshold m t = 5.66. 
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Figure 6: 

Integration domain for calculating inter-earthquake rates in a catalog of duration T. The 
minimum time interval, to, corresponds to the coda duration time of the first event. 



43 



Fig. 7 
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Figure 7: 

Earthquake pair numbers in all zones with rates — 26 x 10~ 20 events/m 2 /s. Earthquake rate 
is taken from a table by Bird et al. (2010) for a magnitude threshold m t = 5.66. Pair number 
is 240. Solid curve is earthquake pair numbers; dashed curve is Poisson approximation ffT5"j) . 
dash-dotted curves are for power-law approximations ffTB"]) . the 6- value is 0.5, 0.65, 0.75, 0.85, 
0.925, and 0.99 from top to bottom. The catalogs time interval is 1976/1/1-2010/11/14. 
Average recurrence time (t) for earthquakes is i — 3272 ± 3596 days. 
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Fig. 8 
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Figure 8: 

Earthquake pair numbers in all zones with rates 26 — 45 x 10~ 20 events/m 2 /s. Pair number 
is 333. Average recurrence time (t) for earthquakes is i — 3507 ± 3242 days. For notation 
see Fig. 
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Fig. 9 
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Figure 9: 

Earthquake pair numbers in all zones with rates > 45 x 10 -20 events/m 2 /s. Pair number is 
264. Average recurrence time (t) for earthquakes is i — 3947 ±3319 days. For notation see 

Fig.m 
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Fig. 10 
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Figure 10: 



Earthquake pair numbers in active continental zones. Pair number is 37. Average recurrence 



time (t) for earthquakes is t = 1304 ± 2578 days. For notation see Fig. 
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Figure 11: 



Correlation between maximum likelihood estimates of fi and 9 parameters in several earth- 
quake catalogs (Kagan et al, 2010). The subscripts i at correlation coefficients (pi) point to 
the table number in Kagan et al. (2010). 
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Fig. 12 
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Figure 12: 



The average aftershock numbers from the PDE catalog in logarithmic time intervals follow- 



ing m > 8.0 GCMT earthquakes during 1977-2003. Line with plus signs shows ra& > 4.9 



aftershocks; dashed lines are theoretical estimates for the first generation aftershocks. Two 



curves display results of earthquake catalog simulation: line with circles shows first genera- 



tion aftershocks, line with squares indicates total aftershock numbers. 
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Fig. 13 

p = 0.83, L = 2.48 + 0.492 (m w -8.25) ' 

a = 0.1 34, e =0.468, n = 160 

max 

L = 2.48 + 0.493 (m w -8.25) + 0.00133 (m w -8.25)' 

a = 0.134,e =0.468 
max 




7.5 



8.5 



Moment Magnitude (m w ) 



Figure 13: 



Plot of log aftershock zone length (L) against moment magnitude (m). Rupture length 
is determined using a 1-day aftershock pattern. Values of the correlation coefficient (p), 
coefficients for linear (dashed line) and quadratic (solid line) regression, standard (a) and 
maximum (e max ) errors, and the total number (n) of aftershock sequences are shown in the 
diagram. 

Circles - thrust mainshocks; 
Stars - normal mainshocks; 
Pluses - strike-slip mainshocks. 
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Aftershocks vs moment, active cont. quakes (m > 7.0, 1977-2010) 
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Figure 14: 



Plot of log aftershock zone length (L) against moment magnitude (m) for earthquakes in 
active continental zones (Kagan et al, 2010). For notation see Fig. 
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